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ABSTRACT 

We present the angular distribution of gamma rays produced by proton-proton interactions in 
parameterized formulae to facilitate calculations in astrophysical environments. The parameteriza- 
tion is derived frorn Monte Carlo simulat ions of the up-to-da te proton-proton interaction model by 
iKamae et al.l (|2005l ) and its extension bv iKamae et al.l (j200 6l). This model includes the logarithmi- 
cally rising inelastic cross section, the diffraction dissociation process and Feynman scaling violation. 
The extension adds two baryon resonance contributions: one representing the A(1232) and the other 
representing multiple resonances around 1600 MeV/c^. We demonstrate the use of the formulae by 
calculating the predicted gamma-ray spectrum for two different cases: the first is a pencil beam of 
protons following a power law and the second is a fanned proton jet with a Gaussian intensity profile 
impinging on the surrounding material. In both cases we find that the predicted gamma-ray spectrum 
to be dependent on the viewing angle. 

Subject headings: cosmic rays — galaxies: jets — gamma rays: theory — ISM: general — neutrinos 
— supernovae: general 



1. INTRODUCTION 

Gamma-ray emission due to decays of neutral pi- 
ons produced in proton-proton {p-p) interactions has 
been predicted from the Galactic ridge, supernova 
remnants (SNRs), active galactic nucleus (AGN) jets, 
and other astronomical sites fiavakawa 1969: Stecken 
19701 iMurthv &: \yolfendaI3 [l986c iSchonfelder, ,200l 
Schlickeiseil 120021: lAharonianl |2004 . A multitude of 
gamma-ray sources a re already known today (see, e.g., 
iHartman et al.l |1999[ ) and with new gamma-ray obser- 
vatories covering GeV to TeV energies many more a re 
expected to be found (lAharonia n et al."2003'. '2004al!bl, 
2005; One 199^: lSchroedter et al. 2005: Schonfelder 2001; 
Weekea i2003i) . The GeV energy window, in particu- 
lar above about 10 GeV, has been poorly explored and 
the GLAST Large Area Telescope (GLAST-LAT)"* is ex- 
pected to provide high-statistics data in this window. 

Interpretation of the observed gamma-ray spectra and 
identification of the involved interactions require not only 
high-quality observational data, but also good knowledge 
of the contributing production mechanisms. In the high- 
energy regime (sub-GeV to multi-TeV energies), the two 
dominant processes for gamma-ray production are p-p 
interactions and subsequent decays of neutral pions, and 
inverse Compton (IC) up-scattering of low-energy pho- 
tons on high-energy electrons. Production of gamma rays 
through pion decay relies on acceleration of cosmic ray 
primaries, protons or heavier nuclei, to high energies. 

The diffuse gamma-ray emission from the Galactic 
ridge was first detected by the OSO-3 and SAS-2 satel- 
lites and later the COS-B and EGRET instrmnents. 
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This emission is interpreted as predominantly due to de- 
cays of neutral pions produced in interactions between 
accelerated protons o r heavier nuclei with the inter- 
stella r medium (ISM; ISteckei" 11973!, IgSp; 'S trong et aD 
197|Jl982, 2000, 2004, 2007; Stephens & Bad hwaijll98'ir 
Dermei1il986atiHunter et al.iil997 ). The gamma-r ay flux 
and spectral shape measured by EGRET (Hunter et al.l 
Il997f) is considered as the key attestation of this inter- 
pretation. There is also significant contribution to the 
diffuse gamma-ray emission from IC scatterings, in par- 
ticular in th e Galactic ridge (Murthv fc Wolfendaldll986t 
IStrong et a l. 2000; Sch5nfcldcr 2001). 

Several SNRs have been detected in TeV energies with 
ground based Air Cherenkov Telescopes (ACTs), in- 
cluding the s hell-type SNRs RX J1713.7-39 46 and RX 
J0852.0-4622 ()Aharonian et al.|[2004al I2005D . Observa- 
tions in the X-ray band from both RX J1713.7-3946 
(iKovama et al.lll997t ISlane et al * 1999; Uchivaina et all 
20031) and RX J0852.0-4622 (.Tsunemi et all 120001 : 
Ivudin et al.ll2005l ) show a smooth, featureless spectrum 
indicating synchrotron X-ray emission due to a popu- 
lation of TeV electrons. The same electron population 
might also produce high-energy gamma rays through 
IC scatterings, but the measured gamma-ray fiuxes and 
spectra in TeV energies do not fully match the pre- 
dicte d ones (see, e.g., the analysis in lUchivama et al.l 
120031 ). and it has been suggested that there might also 
be a significant component du e to hadronic interac- 
tions with the surrounding ISM (iBerezhko fc Volkll2000l: 
Enomoto et"al] 120021 : lAharonianl 120041 : iKatagiri et all 
20051) . Hadronic models fit the very-high-energy gamma- 
ray spectrum asstiming a beam of accelerated protons 
(|Moskalenko et al.l 12007). The highest -energy cosmic 
rays (CRs) escape the forward shock region of the SNR 
almost unidirectionally; the gamma-ray spectrum be- 
comes angular dependent. 

High-energy gamma-ray emission from AGN jets is 
usually explained using leptonic models with electrons 
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accelerated to TeV energies. Support for this comes from 
the observed radio and X-ray spectra which match those 
of synchrotron radiation from high-energy electron pop- 
ulations. The apparent synchronization in the observed 
variability of X-ray and gamma-ray fluxes gives further 
supp o rt for the leptonic niodeling (iOne.1998; Schonfeldeil 
l200lt iSchlickeisei] I200I lAharonianl l200l . There are 
observations of AGN jets where the leptonic scenario 
faces difficulties. For these jets gamma-ray production 
through p -p interactions has been put forward as an al- 
ternative (iMiicke fc Protherod [200lt iMiicke et al']|2003t 
iBottcher fc Reimeill2004l )! 

The giant radio galaxy M87 has recently been ob- 
served i n TeV energies with the H.E.S.S. Cherenkov tele- 
scopes (Bcilickc ct al. 2005). The jet is aligned about 30° 
(|Bickncll & Begclman 1996) relative to the line of sight 
and it has been well studied in radio, optical and X-ray 
wavelengths. The ce ntral object is supposedly a super- 
massive black hole. IStawarz et al.l (j2006) interpret the 
TeV emission as due to Componizat ion of synchro t ron ra - 
diation from a flare in the nucleus. iReimer et al.l (|2004D . 
on the other hand, have suggested the Synchrotron- 
Proton Blazar (SPB) model as the production mecha- 
nism for the TeV gamma-ray emission. The SPB model 
require protons to be accelerated to extremely high en- 
ergies making it less favorable. Another possibility is ac- 
celerated protons leaking out of the jet near the central 
object, interacting with surrounding material and pro- 
ducing gamma-rays through the decay of neutral pions. 

In this paper the angular distribution of gamma rays 
produced in proton-proton interactions is presented in 
parameterized formulae. These are derived from Monte 
Carlo simulatio ns of the up-to-date proton-proton inter- 
action model bv lKamae et al.l ()2006f ). The angular distri- 
bution is given relative to the incident proton direction. 
With this formalism, the gamma-ray spectrum can be 
calculated for any given distribution of protons, includ- 
ing angular dependent ones. 

Parameterization of the angular distributions of other 
stable secondary particles, i.e. electrons, positrons and 
neutrinos, has been deferred due to observational lim- 
itations. When high statistics neutrino data becomes 
available (Halzenll2005h it may be worth while extending 
the parameterization to include the angular distribution 
of neutrinos. 

It is noted that iKoers et al.l ()2006[ ) have taken another 
approach and parameterized the energy and rapidity dis- 
tributions of pions and kaons and from this they are 
able to derive the spectrum and angular distribution of 
gamma rays. They consider proton energies above 1 TeV 
and therefore their model is not suited for studies in the 
GLAST-LAT energy range. In this paper, proton ener- 
gies from the pion production threshold and the reso- 
nance region up to about 10^ GeV are considered. High- 
precision data is expected from the GLAST-LAT in the 
GeV range which makes a parameterization covering this 
range important 

2. PROTON-PROTON INTERACTION MODEL 

For this work the proton-proton interaction model 
used to calcu l ate th e diffuse Galactic ga mma-ray flux 
(I Kamae et al.l l2005l ) and its extension (| Kamae et al.l 
|2006[) is adopted. In an effort to predict the con- 
tribution to the Galactic diffuse emission from 7r° 
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Fig. 1. — Experimental p-p cross sections, as a function of pro- 
ton momentum, and that of adjusted model A: experimental to- 
tal (squares), experimental elastic (triangles), total inelastic (thick 
solid line), non-diffractive (dashed line), diffractive process (dot- 
dashed line), A(1232) (dotted line), and res(1600) (thin solid line). 
The total inelastic is the sum of the four components. The thin 
solid and dot-dot-dashed lines running through the two experimen- 
tal data sets are eye-ball fits to the total and elastic cross sections, 
respectively. The functional f orms are those given by equations (1) 
through (4) in lKamae eTaTI ||2006|) . 



decays, iKamae et al.' (2005D found that past calcu 
lations (iltecker 1970, 1973L [l989l. fStrong et aU lmT 



[Stephens fc Badhwaij[l98ll : Id ermen ll986allbl : lMorilrT99 

left out two important features of the inelastic p-p inter- 
action, the diffraction dissociation process and the Feyn- 
man scaling violation. It was also noted that past calcu- 
lations assumed an energy independent inelastic p-p cross 
section, about 24 mb for Tp ^ 10 GeV, in contradiction 
to recent experimental data where a logarithmic increase 
of the cross section with the incident proton energy is ev- 
ident ('Hagiwa ra et al.ll2002f ). The predicted gamma-ray 
spectrum changed significantly in the GeV energy range 
when the above features were included. The power-law 
index of the gamma-ray spectrum is about 0.05 lower 
in absolute value than that of the incident proton spec- 
trum and the gamma-ray fiux is increased significantly 
comp ared to the reference scaling model ([Kamae ct al] 
I2005D : the increase is proton energy dependent, about 
10-20% in the low-GeV range and about 50% above a 
few 10 GeV. From here on, the model bv I Kamae et al.l 
(|2005f ) is referred to as model A. 

Model A was primarily aimed at studying the diffuse 
emission from the Galactic ridge and thus concerned pro- 
ton kinetic energies well above 1 GeV. It is not accu- 
ra te near the pion p roduction threshold (see Figure 5 
of iKamae et al.l [20051 ). To correct for this and improve 
the ac curacy for lower proton momenta, IKamae et al.l 
()2006[ ) adjusted model A by including contributions from 
two baryon resonance excitation states, A(1232), repre- 
senting the physical A resonance with a mass of 1232 
MeV/c^, and res(1600), representing several resonances 
with masses around 1600 MeV/c^. The term "baryon 
resonance" refers to both nuclcon resonances with iso- 
spin 1/2 and A resonances with iso-spin 3/2. The 
A(1232) decays to a nucleon (proton or neutron) and 
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one pion (7r+, 7r° or tt~; iHagiwara et al.l [2002f) and the 
other resonance, res(1600), is assumed to decay to a nu- 
cleon and two pions. The extension of model A which 
includes the baryon resonances is from here on referred 
to as adjusted model A. 

Other necessary adjustments to model A forced by the 
introducti on of the bary o n reso nance contributions is de- 
scribed in iKamae et al.l ()2006l . section 3) . This includes 
the adjustment of the non-diffractive inelastic p-p cross 
section to accommodate for the resonances while not ex- 
ceeding the total inelastic cross section. The total inelas- 
tic p-p cross section is shown in Figure [1] together with 
the four component cross sections of adjusted model A. 

Due to paucity of experimental data, a-p, p-He and a- 
He have not been included in this work. Near the Earth 
about 7% of the CR flux is a-particles (.Schl ickciscr 200j) 
and the ISM contains about 10% He by number. Both 
the a-particle and the He nucleus can be approximated as 
four individual nucleons. The error from such an approx- 
imation is expected to be less than 10% for high-energy 
gamma rays. Fermi motion of nucleons and multiple nu- 
cleonic interactions in the nucleus affect the pion pro- 
duction near th e threshold and in the resonance region 
(T„ < 3 GeV; ICrawford et~al1 119801 : iMartensson et all 
|2000() . This wiU enhance the pion multiplicity below 100 
MeV. The need for separate treatment of interactions 
such as p-He, a-p, and a-He is acknowledged. 

3. MONTE CARLO SIMULATIONS 

The parameterization of the angular distribution of 
gamma rays is derived from Monte Carlo simulations 
on the ad justed proton-proton interaction model A as 
given in Kamae et al. (2006). Events were generated for 
each of the four components, the non-diffractive inter- 
action, the diffraction dissociation process and the two 
resonance excitation processes in the following way. For 
the non-diffractive interaction in t he high energy r ange 
(Tp > 52.6 GeV) Pythia 6.2 (.Sjostrand et al.ll200ir ) was 
used w ith the option for multipa rton-level scaling vio- 
lation (jSjostrand fc Skandi l2004f ).^ This was comple- 
mented wit h the parameterizatio n of inclusive pion cross 
sections bv lBlattnig et al.l ()2000[ ) in the low energy range 
{Tp < 52.6 GeV). The diffraction dissociation process 
was simulated with a Monte Carlo code by T. Kamae 
(2004, personal communications)® and the resonance ex- 
citation components were simulated with Monte Carlo 
codes by T. Kamae (2005, personal communications). 

For each of the four components mentioned above, 
events were generated for discrete proton kinetic energies 
(0.488 GeV < Tp < 512 TeV) taken from a geometrical 
series 



r„ = 1000 • 2('~22)/2 GeV, 



0, 



,40. 



(1) 



Each proton kinetic energy, Tp, represents a bin cover- 
ing 2-° '^''Tp to 2°-25Tp. The addition of the resonances 
to the model required an increased sampling frequency 
near the pion production threshold and events were also 
generated for Tp = 0.58 GeV and 0.82 GeV. Events were 
not generated for proton energies where the component 
cross section is very small or zero. 



and 



^ See |http://cepa.fnaLg o y/CPD/MCTuningl| 
|http : / / www . phys .utl.edu/~rtield/cdt 

° it is acknowledged tiiat the latest version of Pythia includes 
the diffractive interaction and that the code used here agrees with 
Pythia. The Monte Carlo code is available upon request. 
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Fig. 2. — Experimental invariant cross section at pro- 

duction angle 6cms = 90° and beam energy y/s = 540 GeV for 
(tt"*" -|- 7r~)/2 (filled triangles) and tt" (tilled squar es) measured 
by the UA2 collaboration (p-p collider experiment; [Banner et al.l 
fTMa.[T98l) and beam energy ^ = 53 GeV for (7r+ -|-7r-)/2 (o pen 
circles) measured at the ISR (p-p collider; lAlper et aLlllOTSh to- 
gether with the TT^ invariant cross section calculated from Monte 
Carlo simulations in this work at Tp = 181 TeV (y^ = 582 GeV, 
solid Une) and Tp = 1.41 TeV (v^ = 51.5 GeV, dashed line). 

3.1. Pion Transverse Momentum 

The above described Monte Carlo simulations have 
been verified to agree with experimental data available 
for pions. The inclusive tt" cross section was verified 
to agree with the experimental one and the simulations 
were ensured to reproduce the distributions of pion ki- 
netic energy in the p -p center-of-rnass ( CM) system in 
the resonance region (jKamae et all 120061 ) ■ 

For this work, the 7r° transverse momentum distribu- 
tion was compared with those measured with accelerator 
experiments. Figure [2] shows the invariant 7r° cross sec- 
tion, Ed^a/dp^, at production angle ^cms = 90° for pro- 
ton kinetic energies Tp = 1.41 TeV (y/s = 51.5 GeV) 
and 181 TeV {y/s = 582 GeV) together with experi- 
ment al data for tt^ me asured at the ISR at ^/s = 53 
GeV (jAlper et aLlll975D and bv the UA2 c ollaboration 
at = 540 GeV (|Banner et allfigsl [Ml). One must 
note that the ISR was a p-p collider and that UA2 was 
a p-p collider experiment. The distributions follow the 
expected exponential form for small pt- 

In addition, the energy dependence of the average 
transverse momentum, (pt), calculated from Monte Carlo 
event data was compared with that from ISR experi- 
ments. In accelerator experiments it is difficult to mea- 
sure 7r° directly but one can expect that (p([7r°]) ~ 
(pt[7r±]) (jAlner et al.l Il987t) . (pt[7r°]) was calculated 
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Fig. 3. — Average transverse momentum, {pt), for production 
of pions and all charged particles (tt^, and protons) versus 

the laboratory momentum. Data points are for 7r+ (open trian- 
gles), TT~ (filled triangles) and all charged particles (filled squares) 
from ISR, all charged particles (filled diamonds) from UAl and all 
charged particles (filled circles) from CDF . ISR data are fro m p- 
p collider experiments and are taken from IRossi et all II1975I ) and 
UAl and CDF d ata arc from p-p collider experiments and are taken 
from lAbe eTan (TdSSi) . Lines are for vr" (solid line) and all charged 
particles (dashed line) from Monte Carlo simulations in this work. 

without any fitting and Figure [3] shows the average trans- 
verse momentum as a func tion of the proton momentum 
in the laboratory frame. IRossi et al.l (|1975[ ) estimated 
the error on (pt[7r±]) to be about 10% and {pt[TT°]) cal- 
culated here is within this error margin. 

At very high proton momentum experimental data is 
in general limited to the average transverse momentum 
of charged particles, (p* [charged]), where charged parti- 
cles include charged pions and kaons and protons. Again, 
without fitting, (pt [charged]) was calculated from Monte 
Carlo event data and compared with experimental data 
measured at the ISR and by the UAl and CDF collabo- 
rations.^ As can been seen in Figure [3] the difference is 
on the order of a few percent except at very high ener- 
gies. The discrepancy in average transverse momentum 
at high proton momentum is inherent to Pythia 6.2 which 
has not been fine-tuned to the highest energy data avail- 
able from CDF (T. Sjostrand 2007, personal communi- 
cations). 

3.2. Angular Distribution of Gamma Rays 

The Monte Carlo simulations generate data on mo- 
mentum, p = {px,Py,Pz), and total energy, E, for each 
gamma-ray photon. Since the incident proton direction 
in the simulations is along the z-axis the transverse mo- 
mentum, Pi, is simply 



(2) 



and this was calculated for each simulated event. Events 
were then binned in 2D histograms, 

(3) 

Alog(i?)Apt' ^ ' 

over total energy, E, and transverse momentum, pt, with 
one histogram per proton kinetic energy. Bins widths 
were taken to be Alog(£:) = 0.05 and Apt = 10 MeV/c. 

Note again that ISR was a p-p collider and that UAl and CDF 
are p-p collider experiments. 



Normalization to the proton-proton inela stic cross sec- 
tion cTpp (given by eqs. (1) through (4) in I Kamae et al.l 
[2009) . and per proton-proton interaction gives the differ- 
ential cross section 



A \ogiE)Apt Npp A log{E)Ap 



(4) 



where Npp is the number of proton-proton events simu- 
lated and AiV is the number of gamma rays in a given 
bin. This differential cross section is a representation of 
the angular distribution of gamma rays. 

4. PARAMETERIZATION OF GAMMA-RAY TRANSVERSE 
MOMENTUM DISTRIBUTIONS 

For each proton kinetic energy, Tp, the transverse mo- 
mentum distribution is parameterized as 



A^a 



A\ogiE)Apt 



ptF{pt,x)F\a{pt,x), 



(5) 



where x = log (£'[GeV]), F{pt,x) is the function rep- 
resenting the differential cross section Aa j Apl and 
FwiptlX) is used to approximately enforce the energy- 
momentum conservation. Assuming axial symmetry 
around the p\\ axis, phase space is proportional to 
dpfdp^i^ = 2ptdptdp\\, which gives the extra factor of pt 
in equation ([S]). 

The function Fki(pt,a:) enforcing the energy- 
momentum conservation is taken to be 



1 



exp(VK(pt -Lp)) + r 



F\a{pt:x) 

where W = 7h and 

0.0976 4-0.670 exp(1.81a;) 
_ . -0.793 + exp(0.271(j: + 1) 
^p-\ +0.363(2: -f 1)2) 
2.5 



a; < -1 

-l<x< 0.5 
X > 0.5, 



(6) 



(7) 



with X = log (£;[GeV]). 

In contrast to th e param e teriza tion of inclusive 
cross sections by Kam ae et al.l ()2006[ ). where the non- 
diffractive and the diffraction contribution were treated 
separately, the two are here merged to one contributing 
component. This is well justified in astrophysical con- 
texts. The new component is from here on referred to as 
the non-resonance component. 

The Pt distribution is given by 2D histograms, one his- 
togram per proton kinetic energy, Tp, and component: 
non-resonance, A(1232), and res(1600). Each histogram 
is fitted in slices along pt, i.e. each slice Aa/Ap^ covers 
one bin of Alog(i?). Note that Aa j Apl does not imply 
integrating over \og{E). 

For the non-resonance component Aa j Ap\ is expected 
to follow an exponential form 



Fnript.x) 



ao exp 

ai 



(8) 



Parameter ai gives the shape of the differential cross sec- 
tion and ao gives the absolute normalization. When in- 
tegrating over Pt one should recover the inclusive cross 
section Aa/ A\og{E), i.e. 



Aa 



dpt = 



Aa 



Alog(S)' 



(9) 
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TABLE 1 

Parameters describing transverse momentum distributions 



Parameters Formulae as functions of the proton kinetie energy, y = log (Tp [TeV] ) 

Non-resonance, eg. 11121) 



aio 0.043775 -I- 0.010271 exp(-0.55808i/) 

an 0.8 

ai2 0.34223 + 0.027134j/ - 0.0089229s/2 + 4.9996 x IQ-^y'^ 

ai3 -0.20480 -I- 0.013372y -|- 0.13087 exp {0.0044021(s/ - 11.467)2) 

ai4 ai{x = —0.75) 



A(1232), eq. JlTj 

fcio 18.712 -I- 18.030y + 5.8239y2 -I- 0.62728y^ 

fell 612.61 -I- 404.80y + 67.406y2 

fei2 98.639 -I- 96.741y + 31.597y2 + 3.4567?;^ 

fei3 -208.38 - 183.651/ - 53.283j/2 - 5.0470j/3 

fe2o 0.21977 + 0.064073a' 

fe2i 3.3187 X 10^ + 3463. 4j/ -I- 1.1982 x lO'^y^ + 136.71j/^ 

622 91.410 -I- 91.613y + 30.621y2 + 3.4296j/^ 

623 -521.40 - 529.06y - 178.49j/2 - 19.975?/=* 



res(1600), eq. ((T8J 

cio -1.5013 - 1.12811/ - 0.19813{/2 

cii -33.179 - 22.4961/ - 3.31081/2 

C12 116.44 -I- 122.111/ + 42.5941/2 + 4.9609i/^ 

ci3 -545.77 - 574.801/ - 201.25i/2 - 23.4001/=* 

C20 0.68849 + 0.364381/ -I- 0.047958i/2 

C21 -1.6871 X 10"* - 1.7412 x lO^i/ - 5.9648 x 10^*1/2 - 679.27^=* 

C22 -88.565 - 94.0341/ - 33.014i/2 _ 3.8205i/=* 

C23 1.5141 X 10=^ + 1.5757 x lO^y -I- 544.20i/2 + 62.446i/3 



Thus, ao is taken such that 

/•OO 



ao 



which gives 



Ptexp 



ao 



-— ]dpt = 
ai 



1 Act 
of AlogCE) 



Act 



Alog(£;) 



(10) 



(11) 



and Act/A log(i?) is calculated usi ng the parameteriza - 
tion of the inclusive cross section bv Kamae et al.l(|2006D . 

Parameter ai is a function of both the gamma-ray en- 
ergy, E, and the proton kinetic energy, Tp. It is first 
fitted as a function of a; = log(£'[GeV]) for each simu- 
lated proton kinetic energy. The formula describing ai 
is 

flio exp (-aii(x -I- 012)^) a; < -0.75, 
ai3(x-f 0.75)-|-ai4 a; > -0.75. ^ > 

The parameters an (j = 0, . . . , 4) are then given by func- 
tions of the proton kinetic energy, which are listed in 
Table m 

For the baryon resonance components Act/A^^ will not 
follow the exponential form. Instead, x) is fitted to 

a Gaussian form 

J? I \ v. f (Pt- bi)^ ' 

FA{i232){Pt,x) = 60 exp 



ai(x) 



(13) 



and 



Fi-cs{i600) {Pt,x) = Co exp 



(Pt 



C2 



(14) 



With the requirement that the integral over pt should 
recover the inclusive cross section (eq. [9]) 

60 = 2(6iy^(erf(&i/v/fo^) + 1) + 



+ b2 exp(-62/62))- 



Act 



and 



co = 2(ci0rc^(erf(ci/^) + 1) + 



-|-C2exp(-Ci/C2)) ^ 



Act 



Alog(£;) 



(15) 



(16) 



Again, the parameters bi and Ci (i = 1, 2) are functions 
of both E and Tp and the same procedure is followed for 
them, with 



bi{x) = 6ioexp -bii 



- b 



i2 



1.0 + 5^3(2; -6^2) 



(17) 



and 




Ct{x) = Cjoexp 



for X < 0.5 and bi{x) = for x > Xb and Ci{x) 
X > Xc, with 



Xb^0.8l{y + 3.32) - 0.5 
a;c = 0.82(y-|-3.17) - 0.25, 



(19) 
(20) 



where y = log (Tp [TeV]). These limits of bi and Ci were 
introduced to control artifacts near the kinematical lim- 
its. The parameters 6y and Cy {j — 0, ... ,4) are listed 
in Table [1] as functions of the proton kinetic energy. 

Figure 3] shows the gamma-ray differential cross sec- 
tion Aa/Apt for the non- resonance contribution cal- 
culated using the above described parameterization for 
proton kinetic energies Tp = 64 TeV and 8 TeV and 
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Fig. 4. — Gamma-ray diflferential cross section Aa/Apt for the non-resonance contribution calculated using the parameterization (thick 
solid line) and superimposed with the Monte Carlo simulated cross section (thin histogram). Panels on the left-hand side are for proton 
kinetic energy Tp = 64 TeV and panels on the right-hand side are for Tp = 8 TeV. Rows are for gamma-ray energy E = 1 GeV (top), 100 
GeV (middle), and 1 TeV (bottom). 



gamma-ray energies E — 1 GeV, 100 GeV and 1 TeV. 
The plots show how the average transverse momentum, 
(pi), increases with increasing gamma-ray energy. Fig- 
ure[5]shows Aa/Apt for the two resonance contributions, 
A(1232) and res(1600), calculated at proton kinetic en- 
ergy Tp — 0.82 GeV and gamma-ray energy E — 0.3 
GeV. Superimposed in both figures are the differential 
cross sections from the Monte Carlo simulations. The 
agreement is in general good except near the higher and 
lower kincmatical limits where low statistics in the Monte 
Carlo simulations limits accuracy of the fit. 

5. APPLICATION OF FORMULAE 

The parame terized model of sta ble secondary parti- 
cle spectra by iKamae et al.l (|2006D was used to predict 



differences in the diffuse gamma-ray spectrum from the 
Galactic ridge compared with the scaling models imple- 
mented in GALPROP. The present model finds its ap- 
plication in scenarios where the gamma-ray spectrum is 
expected to be angular dependent, such as AGN jets and 
GRBs, but SNRs may also fall into this category. The 
highest-energy CR escape the forward shock almost uni- 
directionally giving rise to beaming in SNRs. 

To demonstrate the use of the parameterized formulae 
for gamma-ray distributions, the gamma-ray spectrum 
has been calculated for two different cases; the first is a 
pencil beam of protons following a power law of index 
2.0 and the other is a fanned proton jet with a Gaussian 
angular profile impinging on the surrounding matter. 



Parameterization of Gamma-ray Angular Distribution 
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Fig. 5. — Gamma-ray differential cross section Act/ Apt for the 
two resonance contributions, A{1232) and res(1600), calculated us- 
ing the parameterization (thick solid line) and superimposed with 
the Monte Carlo simulated cross section (thin histogram) for pro- 
ton kinetic energy Tp = 0.82 GeV and E = 0.3 GeV. The top panel 
is for for the A(1232) resonance and the bottom panel is for the 
res(1600) resonance. 

5.1. Pencil Beam of Protons 

Consider a beam of protons along the z-axis with no 
spatial extension in the x-y plane, i.e. a pencil beam. 
The energy distribution of protons is assumed to be a 
power law, dN/dE = T~^, with index s — 2.0 and ex- 
tending up to Tp = 512 TeV. The gamma-ray spectrum, 
E^dF/dE, is calculated for three different observation 
angles 6 — 0° (head on), 0°5, and 2° relative to the beam 
axis. The spectra, which are shown in Figure [51 are in- 
tegrated over the annular portion {6, 9 + d9) of width 
d6 = 2'. The absolute normalization is relative to the 
density and distribution of target protons. For compari- 
son the spectrum integrated over the entire phase space 
is also plotted in the figure. As can be seen in the figure, 
the gamma-ray emission is peaked in the very forward 
direction. When the viewing angle is increased, the peak 
of the spectrum is shifted to lower gamma-ray energies 
and the the flux decreases rapidly. Fluctuations in the 
histograms are due to low event statistics. 

5.2. Fanned Proton Jet 

The second example is a fanned proton jet which fea- 
tures a Gaussian intensity profile centered on the jet axis. 
With a FWHM of 3° the opening angle of the jet is about 
10°. The gamma-ray spectrum of the jet is integrated 
over the intensity profile, which is sampled in 0°1 x 0°1 
bins, where each bin is represented by the average of ten 
randomly sampled pencil beams pointing within the bin. 



10-' 1 10' 10' 10' 

|GeV] 

Fig. 6. — Gamma-ray spectra from a pencil beam of protons 
observed from three different angles, 9 = 0° (solid), 0°5 (dashed), 
and 2° (dotted) calculated using the parametric model. The spec- 
tra are integrated over the annular portion (0, 9 + d9) of width 
d8 = 2'. Included is also the spectrum integrated over the en- 
tire phase space (dash-dotted). Histograms are the corresponding 
Monte Carlo spectra. The protons in the beam are assumed to fol- 
low a power law in kinetic energy, Tp, with index 2.0 and extending 
up to Tp = 512 TeV. Fluctuations in the histograms are due to low 
event statistics 




E., |GeV| 

Fig. 7. — Gamma-ray spectra, calculated per solid angle, from a 
proton jet, with a Gaussian intensity profile (FWHM 3°) centered 
on the jet axis, observed from four different angles, 9 = 0° (solid) 
5° (dashed), 10° (dotted), and 20° (dot-dashed). The protons in 
the jet are assumed to follow a power-law distribution with index 
s = 2.0 and extending up to Tp = 512 TeV. 

Protons in the jet are again assumed to follow a power- 
law distribution with index s = 2.0 and extending up to 
Tp = 512 TeV. The gamma-ray spectra, calculated per 
solid angle, observed from four different angles, 9 = 0° 
(head on), 5°, 10°, and 20° are shown in Figure [T] As 
with the pencil beam spectra, the absolute normaliza- 
tion is relative to the density and distribution of target 
protons. 

When the viewing angle is smaller than the opening 
angle of the jet the gamma-ray spectrum features a tail 
extending up to the highest possible gamma-ray energy, 
about 10^ GeV, as can be seen in Figure [T) The tail 
is suppressed for larger viewing angles because of the 
Gaussian intensity profile. At 6* = 5° the tail is about 
four orders of magnitude lower in flux. 

6. CONCLUSIONS AND FUTURE PROSPECTS 
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The angular distribution of gamma-rays produced by 
proton-proton interactions have been presented in pa- 
rameterized formulae. The formulae were derived from 
Monte Carlo simulati ons of the up-t o -date proton-proton 
interaction model by Kam ae et al.l ([2006) and they fa- 
cilitate computation of gamma-ray spectra in cases of 
anisotropic proton distributions. The formulae incorpo- 
rate all important known features of the proton-proton 
interaction up to about Tp — 500 TeV. 

As an example of the application of the formulae, 
gamma-ray spectra was calculated for different viewing 
angles from a pencil beam of protons as well as a pro- 
ton jet with an Gaussian intensity profile. The pencil 
beam example shows very clearly that the gamma-ray 
spectrum changes significantly as the observer is moved 
off the beam. The gamma-ray flux decreases drastically 
and the spectrum gets cut off at lower energy. 

The jet with a Gaussian intensity profile provides a 
more realistic example. The flux does not decrease as 
drastic, but the spectrum changes significantly. For on- 
axis observers the spectrum features a prominent tail 
which is suppressed as the observer is moved off axis. 
The peak of the spectrum is shifted to lower energies as 
the observer is moved. 

Particle acceleration models predict some degree of 

^ http:/ /www. slac.stanford.edu/~niklas/cparamlib 



anisotropy for the highest-energy particles escaping from 
the acceleration site. The parameterized model pre- 
sented here can be used to calculated the anisotropy in 
gamma-ray emission for any given anisotropy in the pro- 
ton distribution. 

The implementation of formulae and parameters given 
in this paper in a C language library will be made avail- 
able by the authors*. The functions implemented in this 
library can be used to calculate both the parameters in 
Table [T] for any given Tp and the differential cross section 
in equation [5] using those parameters for any given set of 
E and pt- 
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